Dealing with Missing Predictor Values

Naturally, we expect that some examples are missing features, let’s try removing those samples and see how much data we’d have left to work with:

train_data_omit_na = na.omit(train_data)
nrow(train_data_omit_na)

## [1] 0

Bummer! It appears every sample has is missing at least one feature. So we can’t simply omit them or we’ll have nothing left to work with! For factor variables with missing values, we can simply create a new other category to fill in for NA values.

other_category = 'other'
for (f in features) {
  if (is.factor(train_data[[f]]) && any(is.na(train_data[[f]]))) {
    levels(train_data[[f]]) = c(levels(train_data[[f]]), other_category)
    train_data[[f]][is.na(train_data[[f]])] = other_category
  }
}

# Try omitting samples with missing features again:
train_data_omit_na = na.omit(train_data)
nrow(train_data_omit_na)

## [1] 1121

This appears to solve most of the problem, but we’re still losing 339 samples, since they have missing numeric feature values. Since they’re a large portion (23% of the dataset), we don’t want to ignore them. Unfortunately, finding suitable replacement values for numeric features is more difficult.

Let’s see which features have missing values:

Filter(function (f) any(is.na(train_data[[f]])), features)

## [1] "LotFrontage" "MasVnrArea"  "GarageYrBlt"

Great, we only need to deal with missing values for 3 numeric features. Let’s replace the NA values of these features with sensible defaults.

For MasVnrArea, the masonry veneer area in square feet, we have a large number of examples with a value of 0, so that should be a decent replacement. Since it doesn’t seem especially important, we might also consider dropping the feature entirely instead of the samples.

train_data$MasVnrArea[is.na(train_data$MasVnrArea)] = 0

Since GarageYrBlt and YearBuilt are highly correlated (0.8235195) and they’re identical 60% of the time, we’ll simply replace missing GarageYrBlt values with YearBuilt values.

train_data$GarageYrBlt[is.na(train_data$GarageYrBlt)] = train_data$YearBuilt[is.na(train_data$GarageYrBlt)]

Coming up with replacement values for LotFrontage (linear feet of street connected to property) is much more difficult. It seems too important of a feature to simply remove. There are also no cases where the value is 0, and so no obvious replacement value.

There are a few variables like LotArea, LotShape and LotConfig that seem like they might be predictive of LotFrontage. Let’s quickly build a model to predict LotFrontage and then use it to fill in our missing values.

lot_frontage_model = lm(LotFrontage ~ log(LotArea) + LotArea:LotConfig + LotShape + LandSlope + YearBuilt + BldgType, data = train_data_omit_na)

We can be confident the lot_frontage_model above is doing a good job predicting LotFrontage because:

• The p-value of the significance of regression test is extremely small.
• The adjusted \(R^2\) is fairly large: 0.583933
• We’re using a small number of predictors without higher degree polynomials so we aren’t especially concerned with overfitting.

If we wanted to be even more confident, we might check the cross validated RMSE, but because we’re only using this model to fill in missing values, we won’t spend too much time on it.

Now we can fillin the missing LotFrontage values:

train_data$LotFrontage[is.na(train_data$LotFrontage)] = predict(lot_frontage_model, newdata = train_data[is.na(train_data$LotFrontage),])

# Check if there are any remaining NA values in the training data:
any(is.na(train_data))

## [1] FALSE

Success! We’ve managed to replace all of the missing values without losing any training samples.

Creating Additional Factor Variables

When we loaded the dataset, R decided which predictors were numeric and which were factors. Some predictors that R considers numeric, might work better as factor variables. For example, the price difference between a 0 and a 2 car garage is likely different than the price difference between a 2 and a 4 car garage.

for (f in features) {
  # We'll add any numeric predictor with less than 12 unique values as a factor variable
  if (!is.factor(train_data[[f]]) && length(unique(train_data[[f]])) <= 12) {
    train_data[[paste(f, 'Fctr', sep = "")]] = as.factor(train_data[[f]])
  }
}

# Update `features`:
features = setdiff(colnames(train_data), c("Id", "SalePrice"))

Great, we’ve added 14 new factor variables for a total of 93 features.

Removing Categories with Few Samples

// TODO: explain this

category_cutoff = n * 0.05
for (f in features) {
  if (is.factor(train_data[[f]])) {
    for (category in levels(train_data[[f]])) {
      samples_in_cat = sum(na.omit(train_data[[f]]) == category)
      if (samples_in_cat < category_cutoff) {
        levels(train_data[[f]]) = c(levels(train_data[[f]]), other_category)
        train_data[[f]][train_data[[f]] == category] = other_category
      }
    }
  }
}

feature_cutoff = n * 0.95
for (f in features) {
  if (is.factor(train_data[[f]])) {
    if (length(levels(train_data[[f]])) == 1) {
      train_data = train_data[,!(names(train_data) == f)]
    } else {
      for (category in levels(train_data[[f]])) {
        samples_in_cat = sum(na.omit(train_data[[f]]) == category)
        if (samples_in_cat > feature_cutoff) {
          train_data = train_data[,!(names(train_data) == f)]
        }
      }
    }
  }
}

Baseline Models

# Start with a simple, naive model using all predictors added together 
additive_large = lm(SalePrice ~ ., data = train_data)
additive_large_eval = evaluate_model("Additive Large", additive_large)

## [1] "[cache] loaded train model from  ~/git/STAT420/final-project/cache/1dc7dcfe2294b58c7db2d8267acde55d53154dc1.rds"
## [1] "[cache] loaded train model from  ~/git/STAT420/final-project/cache/e48f7dab9a7fa924a851ed73636655b162e8129d.rds"
## [1] "[cache] loaded train model from  ~/git/STAT420/final-project/cache/acbaf62b1e90944ba172fd50696957ad44848541.rds"
## [1] "[cache] loaded train model from  ~/git/STAT420/final-project/cache/c77232ccfd0247067ea5946281d8b361174baa27.rds"
## [1] "[cache] loaded train model from  ~/git/STAT420/final-project/cache/c77232ccfd0247067ea5946281d8b361174baa27.rds"

kable(additive_large_eval)

# Variables were selected based on the results of the Data Exploration section above
additive_small = lm(SalePrice ~ 
                      LotArea + X1stFlrSF + X2ndFlrSF + MasVnrArea + BedroomAbvGr + OverallQual + 
                      OverallCond  + ExterQual + KitchenQual + BsmtQual + BsmtExposure + ScreenPorch + 
                      Fireplaces + MSZoning + YearBuilt, data = train_data)
additive_small_eval = evaluate_model("Additive Small", additive_small)

## [1] "[cache] loaded train model from  ~/git/STAT420/final-project/cache/1dc7dcfe2294b58c7db2d8267acde55d53154dc1.rds"
## [1] "[cache] loaded train model from  ~/git/STAT420/final-project/cache/e48f7dab9a7fa924a851ed73636655b162e8129d.rds"
## [1] "[cache] loaded train model from  ~/git/STAT420/final-project/cache/acbaf62b1e90944ba172fd50696957ad44848541.rds"
## [1] "[cache] loaded train model from  ~/git/STAT420/final-project/cache/c77232ccfd0247067ea5946281d8b361174baa27.rds"
## [1] "[cache] loaded train model from  ~/git/STAT420/final-project/cache/c77232ccfd0247067ea5946281d8b361174baa27.rds"

kable(additive_small_eval)

We can do much better on our target metric (CV.Log.RMSE) with the smaller model. Notice that the small additive model has far fewer parameters than the large one (25 vs 172) It seems that with so many predictors, it’s very easy to overfit to the training data.

Parameter Search using BIC

The small baseline model also has a better BIC. So it looks like BIC might be a good metric to use to avoid overfitting. Next, we’ll use BIC to do a backwards parameter search starting from additive_large model.

additive_bic_backward = cached_step(additive_large, k = log(n), direction = 'backward', trace = 0)

## [1] "[cache] loaded step model from  ~/git/STAT420/final-project/cache/923e6ff8a018aa605fe3e56572395ef5434e8a71.rds"

additive_bic_backward_eval = evaluate_model("Additive Backward BIC", additive_bic_backward)

## [1] "[cache] loaded train model from  ~/git/STAT420/final-project/cache/1dc7dcfe2294b58c7db2d8267acde55d53154dc1.rds"
## [1] "[cache] loaded train model from  ~/git/STAT420/final-project/cache/e48f7dab9a7fa924a851ed73636655b162e8129d.rds"
## [1] "[cache] loaded train model from  ~/git/STAT420/final-project/cache/acbaf62b1e90944ba172fd50696957ad44848541.rds"
## [1] "[cache] loaded train model from  ~/git/STAT420/final-project/cache/c77232ccfd0247067ea5946281d8b361174baa27.rds"
## [1] "[cache] loaded train model from  ~/git/STAT420/final-project/cache/c77232ccfd0247067ea5946281d8b361174baa27.rds"

kable(additive_bic_backward_eval)

Transformations

Doing a parameter search using backwards BIC managed improve our results, but it still fails the Shapiro-Wilks and Breusch-Pagan tests. This suggests we might be able to achieve better results by applying some transformations to our the response and/or predictors.

Based on the graphs in the Data Exploration section above, it appears that many predictors have an exponential or polynomial relationship with SalePrice.